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Abstract 

Extremal dynamics is the mechanism that drives the Bak- 
Sneppen model into a (self-organized) critical state, marked by 
a singular stationary probability density p(x). With the aim of 
understanding of this phenomenon, we study the BS model and 
several variants via mean-field theory and simulation. In all cases, 
we find that p(x) is singular at one or more points, as a con- 
sequence of extremal dynamics. Furthermore we show that the 
extremal barrier cc, always belongs to the 'prohibited' interval, in 
which p(x) = 0. Our simulations indicate that the Bak-Sneppen 
universality class is robust with regard to changes in the updating 
rule: we find the same value for the exponent tt for all variants. 
Mean-field theory, which furnishes an exact description for the 
model on a complete graph, reproduces the character of the prob- 
ability distribution found in simulations. For the modified pro- 
cesses mean-field theory takes the form of a functional equation 
for p(x). 
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I. INTRODUCTION 



The Bak-Sneppen (BS) model was proposed as a possible explanation of 
mass extinctions observed in the fossil record [1] , and was recently adapted to 
model experimental data on bacterial populations [2,3]. Independent of its bi- 
ological interpretation, the model has atracted much attention as a prototype 
of self-organized criticality (SOC) under extremal dynamics [4,5]. The model 
has been studied through various approaches, including simulation [6-8], theo- 
retical analysis [9-11], probabilistic analysis (run time statistics) [12,13], renor- 
malization group [14,15], field theory [16] and mean-field theory [4,5,17]. Some 
variants have been proposed, for example the anisotropic BS model [18,19]. In 
this paper, we study the consequences of extremal dynamics, using mean-field 
theory and simulation. With this aim, we propose variants of the model and 
analyse how varying the updating rule affects the stationary probability den- 
sity and the critical behaviour. 

In the evolutionary interpretation of the BS model, each site % represents a 
"niche" occupied by a single species, and bears a real- valued variable Xi repre- 
senting the "fitness" of this species. (In the present context "fitness" denotes 
a propensity to resist extinction: if x^ > Xj then species j goes extinct before i, 
so that Xi may be termed a "barrier to extinction" .) At each step the site with 
the smallest x^, and its nearest neighbors, are replaced with randomly chosen 
values. The replacement of the neighboring variables with new random values 
may be interpreted as a sudden unpredictable change in fitness when a nearby 
niche (which might have borne a predator, or a food source for the species 
in question), is suddenly colonized by a new species. Selection, at each step, 
of the global minimum of the {x^ ("extremal dynamics") represents a highly 
nonlocal process, and would appear to require an external agent with complete 
information regarding the state of the system at each moment. Applications 
of the model in evolution studies are reviewed in [20]. 

In the original updating rule, neighbors are randomly affected by the ex- 
tinction of an interacting species. There is no a priori reason to expect that 
evolution should obey this specific rule on a specific lattice. Thus, we ask: 
what happens if the extinction of the least adapted species favors (or pre- 
vents) the extinction of other species? Moreover, what is the signature of 
extremal dynamics that can, in principle, allow us to recognize it in the real 
world? 

Due to extremal dynamics, the BS model exhibits scale invariance in the 
stationary state, in which several quantites display power-law behaviour [1]. 
Simulations show that the stationary distribution of barriers follows a step 
function, being zero (in the infinite-size limit) for x < x* ~ 0.66702(8) [6]. 
Relaxing the extremal condition leads to a smooth probability density and 
loss of scale invariance [4,5]. 

A striking feature of the BS model is that a simple updating rule leads to a 
singular stationary probability density p(x). An intriguing question therefore 
arises, as to how changes in the updating rule affect this density, an issue that 
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has not, to our knowledge, been investigated previously. In this work we exam- 
ine the consequences of rules in which one or more sites are updated according 
to x — > x' = f(x), instead of being replaced by a random number. (Here 
/ maps the interval [0,1] to itself.) We find that this can provoke dramatic 
changes in the stationary probability density. In the extremal dynamics limit, 
the variant models belong to the same universality class as the original. We 
find that the hallmarks of extremal dynamics are that i) the stationary prob- 
ability density is singular, and ii) with probability one, the extremal Xi (i.e., 
the next variable to be updated) belongs to the 'prohibited' region in which 
p(x) = 0. Using a two-site mean-field approximation, we also find evidence 
that nontrivial correlations are restricted to the prohibited region. 

The balance of this paper is organized as follows. Section II introduces the 
models, which are then analyzed using mean-field like approaches in Sec. III. 
In Sec. IV we present simulation results, and summarize our conclusions in 
Sec. V. 

II. MODELS 

The Bak-Sneppen [1] model is a discrete-time Markov process on a d- 
dimensional lattice of L d sites, with periodic boundaries. At each site we 
define a real- valued variable Xj(0). Initially, these variables are independently 
assigned random values uniform on [0,1). At time 1, the site m bearing the 
minimum of all the numbers {xj(0)} is identified, and it, along with its 2d near- 
est neighbors, are given new random values, again drawn independently from 
the interval [0,1). (In the one dimensional case considered here this amounts 
to: x m (l) = 7], x m+ i(l) = rj', and x m -i(l) = rj", where rj, r/, and r/" are inde- 
pendent and uniformly distributed on [0,1); for \j — m\ > 1, Xj(l) = Xj(0).) 
At step 2 this process is repeated, with m representing the site with the global 
minimum of the variables {^(l)}, and so on. In the random neighbor version 
of the model, the process is realized on a complete graph (all sites are con- 
sidered neighbors); two randomly selected sites are updated in addition to the 
minimum m. 

We now define three modified Bak-Sneppen models, that differ from the 
original only in the way that the barriers X{ evolve. In one, the site M bearing 
the maximum of the {xj} is replaced with a random number r] and the two 
nearest neighbors are replaced with their own square roots: Xm{1 + 1) — V and 
XM±i(t + 1) = \JxM±i(t). We shall refer to this as the 'radical' variant. In 
the second variant, the site with the maximum value is replaced by its own 
value squared, while its two nearest neighbors receive random numbers 77 and 
rj': XM(t + 1) = [xM(t)] 2 , XM+i(t + 1) = ?7, and x M -i(t + 1) = rj'. This will 
be called the 'centered square' version. Finally, we define a 'peripheral square' 
variant, in which one of the nearest neighbors of M is squared, while M and 
its other neighbor are replaced with random numbers: XM'(t + 1) = [ijf'Wf) 
x M (t + 1) = 77 and x M "(t + 1) = 77'. (Here M' = M + a and M" = M - a, 
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where a is a random variable that assumes values of +1 and -1 with equal 
probability.) 

The motivation for studying these variants is twofold. First, the Bak- 
Sneppen model is notable for exhibiting a singular stationary probability den- 
sity, and it is of interest to examine the effect of changes in the dynamical 
rule on this density and on the critical behavior. If we introduce a determin- 
istic function f(x) as part of the updating rule, it is desirable that / map 
the interval [0,1] to itself, making functions of the form f(x) = x a a natural 
choice. In this context we note that the variants feature what may be called 
'migration', that is, the systematic movement of certain variables Xi within 
the interval [0, 1). In the radical variant the migration is from the populated 
region (smaller x) to the 'excluded' region (larger x), whereas in the square 
versions migration occurs in the opposite direction. 

Secondly, the variants admit interpretation as evolutionary processes. (In 
the modified models, we have for convenience defined the site with the maxi- 
mum variable as the most vulnerable, so that small Xi now corresponds to high 
fitness.) In the radical variant, replacement of the least-fit species provokes a 
reduction in the fitness of its neighbors, without leading to their immediate 
extinction. Thus, some memory of the fitness of the neighboring species is 
retained. The radical variant therefore seems a plausible modification of the 
original model, in the biological context. In the peripheral square variant the 
extinction of the least-fit species provokes extinction of one neighbor, and in- 
creased fitness of the other. Finally, in the centered-square variant, xm x 2 M 
represents an increase in the fitness of the least viable element of the system, 
while its neighbors go extinct. 



We develop mean-field approximations for the original and modified mod- 
els, along the lines of Refs. [4] and [5]. To begin, we relax the extremal condition 
introducing a flipping rate of Ye~^ Xl at site i, where is a characteristic time, 
irrelevant to stationary properties, and which we set equal to one (r = 1). Call 
this regularized system the 'finite-temperature' model. The extremal dynamics 
of the original model is recovered in the zero-temperature limit, (3 — > oo. 

Consider the probability density p(x). In the finite-temperature version of 
the original model, p(x) satisfies 



where p(x, y, t) is the joint density for a pair of nearest-neighbor sites and 
p(y, t) is the one-site marginal density. Invoking the mean-field factorization 
p(x,y) = p(x)p(y) (we suppress the time argument from here on), we find: 



III. MEAN-FIELD THEORY 



A. Original Model 



dp(x, t) 
dt 



-e-^p(x, t) - 2 f 1 e-top(x, y, t)dy + 3 f e~^p{y, t)dy , (1) 

JO JO 
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^ = _[ e -/fa + 2/(^)]p(x) + 3/()9) , (2) 

where 

I((3) = [ 1 e-^p(y)dy , (3) 
Jo 

represents the overall flipping rate. Eq. (2) is a nonlinear equation for p, in 
which the density at x is coupled to p at all other arguments via I(/3). 
In the stationary state we have 

p(x) = 3J a . (4) 

Multiplying by e~^ x and integrating over the range of x, we find 

I(/3) = ( e W 3 - l)/[2e p (\ - e~^ 3 )}, (5) 

and thus 

3 i _ e -2/3/3 

Pst(x) = - x _ e _ 2m + e _ px ( em _ x) • (6) 

This solution is plotted for various (3 values in Fig. 1. 

In the limit /5 — > oo the solution becomes a step function: 

p st (x) = ^Q(x - 1/3)6(1 - x) . (7) 

Thus the mean-field approach correctly predicts a step-function singularity for 
the probability density, although it places the critical barrier at a;* = 1/3, 
whereas it actually falls at 0.66702(8) [6]. On the other hand, the slightest 
relaxation of the extremal condition destroys the singularity [5], since p(x) is 
a smooth function for (3 < oo. The rate of convergence to the step-function 
is generally exponential with (3, away from the discontinuity. The curves for 
various (3 values exhibit an approximate crossing near x = 1/3. The derivative 
at this point however diverges only linearly with (3: (dp st /dx) x= i/ 3 ~ 3/3/8 for 
large (3. 

Using Eq. (2), we find, in the limit of large (3, that the relaxation time for 
a small disturbance from the stationary solution grows ~ e' 3 / 3 . (By 'small' we 
mean I(f3) ~ e^/ 3 /2.) 

The following observation will prove useful in the discussion of the modified 
models. If we assume that, in the limit (3 — > oo, p s t{x) is identically zero for 
x < x*, and that p st > C > on some interval [x*,a] (in other words, the 
density suffers a jump discontinuity at x*), then / = Jq e~ /3y p(y)dy ~ e -P x * 
and so e~^ x jl ~ e -P( x ~ x *) _ > o for x > x*. Then Eq. (4) reduces to the step- 
function expression, Eq. (7). Note however that lim^oo e px * I{(3) = 1/2 not 
3/2/3, as would be found by naively inserting the limiting density, Eq.(7), in 
Eq.(3). This means that the dominant contribution to I is due to the interval 
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[0, x*], even in the limit (5 — ■> oo, which is readily seen if we write Eq.(6), for 
large f3, as 

p st (x) ~ |[6(a; - x*) + e(x* - x)e^ x - x ^\ . (8) 

In the limit /3 — > oo, sites with rr < a;* constitute a set of probability zero, but 
the site m selected for extinction belongs to this set with probability one. This 
is a singular property of the Bak-Sneppen model in the infinite-size limit, as 
discussed in the next subsection. (The infinite-size limit is implicit in mean- 
field theory.) 

Although in the modified models we are unable to find an analytical solu- 
tion for finite j3, it is possible to integrate the mean-field equation numerically. 
Due to the factor e^ x , for large f3, a very small time step would be needed 
to avoid instability in the usual integration methods (e.g., Euler or Runge- 
Kutta). We circumvent this difficulty using a partial integration method [21]. 
To apply this method to the MF equation for the original model, we write Eq. 
(2) in the form 

^l = -f(t)p(x,t)+9(t) , (9) 

where /(*) = e~ Px + 21 (t) and g(t) = 31 (t). The formal solution is 

rt' 



git 1 )) • (io) 



p(x,t) =exp[- / duf(u)]\p(x,0)+ [ dt'exp [ dt"f(t") 
Jo [ Jo \_Jo 

For a small time interval At, we find 

p(x, At) ~ e -/(°) A * S^p(x, 0) + g(0) dt'e f W j 

= e-™" P (x,0) + ^(l-e-™") . (11) 

This relation can be iterated to find the evolution of p(x, t) from a given initial 
distribution, which converges to the stationary density. 



B. BS model on a finite complete graph 

Mean-field theory is exact for the "random neighbor" model, which may 
also be thought of as the BS model on a complete graph, i.e., one in which all 
pairs of sites are neighbors. (When x m is updated, two of these neighbors are 
chosen at random for updating as well.) In this subsection we analyze the BS 
model with extremal dynamics on a complete graph of N sites. 

Since sites are assigned independent random numbers, the Xi are indepen- 
dent, identically distributed random variables drawn from the density p(x,t). 
Define the distribution function P(x,t) = J X p(y,t)dy. The probability that 
the next site to be updated, x m , lies between zero and x is: 
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Prob[x m < x] = 1 - [1 - P(x)f , (12) 

i.e., one less the probability that the minimum is larger than x. The probability 
that a randomly chosen neighbor has Xi < x is simply P(x), and the probability 
that one of the updated sites receives a number < x is x. At each step, 
therefore, the expected change in the number of sites with Xi < x is 3x — 
2P(x) — Prob[x m < x], which implies 



dP(x,t) 
dt 



-{l-[l-P(x)] N }-2P(x) + 3x . (13) 



(Here we have taken the time unit to represent N updates.) 
Letting Q = 1 — P, we have in the stationary state, 

Q N + 2Q -3{l-x) = 0. (14) 

(Note that Q(0) = 1, Q(l) = 0, and dQ/dx < 0.) Numerical solution (Fig. 
2) shows that for large N, Pjv(x) approaches a singular function that is zero 
for x < 1/3, while for x > 1/3, P(x) = 3(x — l/3)/2. It is straightforward to 
show that for fixed x and iV — > oo, 

f 1 — (1 — 3^)^ , x < 1/3 
~ < (15) 

l^+l[l(i-x)r, x>i/3. 

For i = 1/3 we have InP ~ — In iV + lnln(A^/2) (plus terms of lower order in 
N) as N — > oo. Of interest is the exponential convergence of P to its limiting 
form for x > 1/3, compared with algebraic convergence for x < 1/3. Note also 
that Prob[x m < 1/3] ~ 1 — 2/N, so that the minimum indeed belongs to the 
excluded region with probability one, when iV — > oo. 



C. Pair approximation 

The analysis of the finite-temperature model is readily extended to the 
pair level, in which one studies the evolution of the two-site joint probability 
density p(x,y). Our starting point is the following exact relation, obtained 
using the same reasoning that led to Eq. (1): 



dp(x, y) 
dt 



- (e-^ + e"^) p(x, y) - jf * e^ u \p(x, y, u) + p(u, x, y)\ du 



+ J J (e- pu + e- pv )p(u,v)dudv 

+ / / [p(x, u, v) + p(v, u, y)} e'^dudv . (16) 
Jo Jo 



Now invoking the pair factorization, 



p(x,y)p(y,u) 

p(x, y, u) ~ — , (17) 

PKV) 
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Eq. (16) becomes 

dp(x, y) 
dt 



o 



e -px + e -py + x(x) + K(y)] p(x, y) + 21 
\p(x, u) + p(y, u)] K{u)du , 



with 



K(x) = J(x)/p(x), (19) 
J(x) = f 1 p( x ,u)e- (Su du (20) 

JO 



and 



/= f p( x )e- /3x dx= f 1 J(x)dx . (21) 

JO JO 

To find the stationary solution numerically, we note that 

_ 2/ + Jo 1 [p(ar, u) + u)] K(u)du 
P ^ y) ~ e -P* + e-0y + K(x) + K(y) ' 1 j 

Starting from an arbitrary normalized density po (for example, uniform on 
[0, 1] x [0, 1]), we generate p\ by evaluating the r.h.s. of Eq. (22) using p and 
normalizing the resulting expression. This procedure is then iterated until it 
converges to the stationary density. We find that for large /3, the stationary 
marginal density approaches the step-function solution 

p ^ = rh^ e{x ~ x * )e{1 ~ x) (23) 

with x* ~ 0.47186, a considerable improvement over the site approximation. In 
this limit, the joint distribution is the product of two identical one-site distri- 
butions. Once again, the portion of the unit square that has probability zero is 
in fact responsable for all transitions. In the region D = {(x, y)\0 < x,y < x*}, 
the two variables are correlated, as shown by the nonzero correlation coeffi- 



cient p = cov(x, y)/\J var(x)var(y). In the pair approximation, p ~ 0.327 for 
(x,y) G D, as (3 — > oo. 

D. Radical Variant 

In the radical model, the probability density p(x) = px(x) satisfies: 
dpx(x) 



dt 



= -e^p x (x)-2 [ 1 e^p x (x,y)dy+ P e^p x (y)dy 
JO jo 

+ 2p xl/2 (x) f 1 e^p x (y)dy , (24) 

JO 
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where Px l / 2 i x ) = 2xpx(x 2 ). With the definition I (ft) = Jq e l3y p(y)dy and the 
mean-field assumption p(x,y) = p(x)p(y), Eq. (24) reduces to 



dp(xj h 



dt 



e px p(x) + I((3)[-2p(x) + 1 + 4xp(x 2 )] . (25) 



Given the step-function form of the stationary density for the original 
model, it is reasonable to expect that in this case as well, p(x) will have a 
jump discontinuity for f3 — > oo, and be zero for x > x*. We can get some 
insight into the nature of the density as (3 — > oo by first observing that if 
p(x) < C < oo as x — > 1, then I ~ Ce 13 / ' (3 for large (3, and therefore e~^I — > 
as (3 — > oo. This in fact holds unless has a 5-like contribution at x — 1, 
which is not expected since it is precisely the largest values of x that are 
removed in the dynamics. Write the stationary solution to Eq. (25) as 

. . 7fl + Axp(x 2 )} . . 

pO»0 = 1 1 + 2 j h , (26) 

with / = e~^ x I. This gives p(l) = in the limit (3 — > oo. Similarly, supposing 
that xp(x 2 ) — > as rr — > 0, we have p(0) = 1/2 in this limit. A similar line of 
reasoning can be used to show that dpjdx\ x= \ = in the (3 — > oo limit. 

The preceding discussion suggests that in the limit (3 — > oo, p s t(x) is «c?en- 
tically zero over some finite range [x*, 1]. Assuming this to be so, we have 
I ~ e^ x and in the limit f3 — > oo the stationary solution is given by the 
functional equation: 

2p(x) - 4xj9(a; 2 ) = 1 , (27) 

for < x < x*. Writing this as p(x) = 1/2 + 2xp(x 2 ), we can iterate to find 
p{x 2 ) = 1/2 + 2x 2 p(x 4 ), p(x 4 ) = 1/2 + 2x 4 p(x 8 ), p(x 8 ) = 1/2 + 2x 8 p(x 16 ), etc., 
which suggests the solution 

= -+x + 2x 3 + 4a; 7 + 8x 15 + ... = ^2 i " 1 a; 2 ^ 1 . (28) 

Substituting this 'lacunary series' in Eq.(27), one readily verifies that pi(x) 
is a solution. Similarly, rewriting Eq.(27) as p(x 2 ) = —l/4:X + p(x)/2x, one 
finds p(x) = -l/Ax 1 ' 2 + p(x^ 2 )/2x^ 2 , p(x 1 ' 2 ) = -l/Ax 1 / 4 + p(x 1 ^)/2x 1 / 4 , 
etc., leading to a second solution: 

1111 00 

V J X ) = = _ V 2- {i+1) x- (1 - 2 ~ t) 

P2 ^ X) ~ Ax 1 / 2 8arV 4 \Q X V 8 32x 15 / 16 "' ^ 

(29) 

We now search a solution of the form p(x) = Api(x) + Bp2(x); substituting 
in Eq. (27), yields the condition A + B = 1. This linear combination must 
however be normalizable. The relevant integrals are: 
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rx* i oo i 

/ Pl (x)dx = - y>* 2 " = -{x* +x* 2 + x* 4 + x* s + ...) , (30) 
Jo 2 n=0 2 



X* 1 OO 1 

p 2 ( x )^ = -W x* 2 ~ n = ~(1 + ** 1/2 + ^* 1/4 + ^* 1/8 + ...) • (31) 

2 n=l 2 



Since < x* < 1, the first integral converges while the second diverges, so that 
P2{x) is not normalizable, implying B — 0. A normalized, positive solution 
is p(x) = pi(x) for x < x* = 0.793189 and p(x) = for x > x*. (x* is 
determined by normalization.) This solution, for infinite /3, is plotted in Fig. 
3. Numerical integration of Eq. (25) through the method outlined in Sec. 
IIIA yields results consistent with this solution, as may again be seen in the 
figure. Finally, simulation of the random-neighbor version of the model yields 
a stationary distribution consistent with this expression (see Fig. 4). 

We will now analyze the radical model on a finite complete graph and show 
that its probability density (equation 28) can be derived via a different path. 
Define the distribution function Q(x,t) = p(y,t)dy. By the same reasoning 
developed in section III B, the probability that x m , the next site to be updated, 
lies between x and 1 is 



Prob[x m > x] = 1 - [1 - Q{x)f . (32) 

The probability that a randomly chosen neighbor has X{ > x is Q(x), while the 
probability that a neighbor receives a barrier > x is Q(x 2 ). The probability 
that x m receives a new value between x and 1 is simply 1 — x. The expected 
change in the number of sites with Xi > x is (1 — x) + 2Q(x 2 ) — 2Q(x) — {1 — 
[1 - Q(x)] N }, so that 

= (1 - x) + 2Q(x 2 ) - 2Q(x) - {l - [1 - Q(x)] N } . (33) 

In the sationary state, the definition P = 1 — Q leads to 

x + 2P(x 2 ) - 2P(x) - P(x) N = . (34) 

(Note that P(0) = and P(l) = 1). Since < P(x) < 1 for x < x*, 
lim^v^oo P(x) N = and the probability density in an infinite system obeys the 
functional equation x + 2P(x 2 ) — 2P(x) = 0. Using the iterative approach, we 
find two solutions, one divergent (which is rejected), and the finite solution 

1 oo 
Z n=0 

i.e., the integral of pi(x), Eq. (28). Normalization then demands that P = 1 
for x > x*. 

Numerical solution of equation (33) (see Fig. 5) shows that, for large N, 
Pn(x) approaches the singular function described above. For fixed x > x* 1 ^ 2 , 
and iV — > oo, we find P(x) ~ 1 — x 1 ^ . 
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E. Centered Square Variant 



The probability density in the finite-temperature version of the centered 
square variant obeys: 



= -e px p x {x) - 2 / e py p x {x, y)dy + 2 / e Py p x (y)dy + —^p x (^/x~) ■ 
Jo Jo l^/x 

(36) 



Mean-field factorization leads to 

^ = -e^p(x) - 21 (P) \p(x) - 1] + ^Lp( yfx) , (37) 
with I(/3) as defined in Sec. HID. In the stationary state, we find 



2 + T±:p(y/x)eW/I 

^ - iZ'/i ■ (38) 

The hypothesis that in the limit (3 — > oo , the stationary density p(x) = for 
x > x* implies I(/3) ~ e" x . Therefore 

,. el3x I if x<x* fon . 

hm — = { ., „ (39) 

/3^oo I [ OO if X > X V ' 

and 

hm — — = <^ *, +2 (40) 

/3^oo I [ OO if X > X* V ' 

This, combined with Eq. (38), implies that p(x) = 1 for x G [0,x* 2 ], and 
that 



p(x) = lim 

/3-+oo 



1 + 



1 



4v^ 



—p(Vx) 



for x G [x* 2 , x*] (41) 



The iterative process used in the preceding section is not useful here as it leads 
to pathological solutions, due to the divergent ratio e 13 ^ J I in this interval. 
Observe that if, for x G [x*, x* 1//2 ], p(x) = g{x 2 )I /e@ x , as (3 — > oo, with g(x 2 ) 
finite, then p(x) = 1 + g(x)/A^Jx for x G [x* 2 ,a;*]. (Note that this means 
that p(x) — > as (3 — > oo for x in the interval [rr*, a;* 1 ^ 2 ], consistent with 
the hypothesis that p(x) — > for x > x* .) The function g(x) is however 
yet to be determined. Attempting the simplest solution, g(x) = constant, 
we find a surprisingly reasonable result, as shown in Fig. 6. Next, allowing 
g(x) = ax + b, with a and b constant, yields excellent agreement with the 
numerical integration of Eq. (37) as also shown in Fig. 6. We do not have an 
argument why g{x) should take this form. 
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The threshold x* can be determined in a simple way. Since the fraction 
of barriers in the interval [x*, 1] is constant in the stationary state, the mean 
number of barriers removed from this interval at each time step must equal 
the mean number inserted. The probability that the maximum x M lies in 
[x*, 1] is 1, while its random neighbors are certainly below x*. Each updated 
neighbor has a probability (1 — x*) of receiving a barrier in the interval [x*, 1], 
while the maximum remains in [x*, 1] with probability p±, the probability that 
% £ [x* 1 / 2 , 1]. Thus, we have 

1 =2{l-x*)+p 1 . (42) 

This reasoning can be repeated for the intervals [x* 1//2 ,l], [x* 1//4 , 1], 
leading to 

Pl = 2(l-x* 1/2 )+p 2 , (43) 

p 2 = 2(l-x* 1/i )+p 3 ... (44) 

where p n is the probability that xm £ [x* 1 ^ , 1]. Substituting this result in 
equation (42), we find 

l = 2(l-a;*) + 2(l-a;* 1/2 ) + 2(l-a;* 1/4 ) + 2(l-a;* 1/8 ) + ... , (45) 

which provides x* = 0.761072. Finally, we note that normalization implies 
f**2 g(x)/4y/x dx = 1 — x* = 0.238928, providing a constraint on the function 
g(x). 



F. Peripheral Square Variant 

We now apply the mean- field analysis to the peripheral square variant. In 
this case, the probability density satisfies: 



<1PxU) = -e^ Px (x) -2 f 1 e^ Px (x,y)dy + 2 \\^p x {y)dy 

Jo Jo 

+ p X 2(x) f 1 e Pv p x {y)dy , (46) 
Jo 



dt 



where px 2 (x) = px(x 1 / 2 )/2x 1 / 2 . Under the mean-field factorization this re- 
duces to 

= -e^p(x) + [-2p(x) + 2 + ^t^ 1/2 )] , (47) 

with I(/3) as given in Sec. IIIC. The hypothesis that, for f3 — > oo, the stationary 
density p(x) = for x > x* implies I ((3) ~ e^ x *, leading to the functional 
equation: 
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for x < x*. 

As in the centered-square variant, the iterative method does not yield a 
useful solution, and we pursue a different approach. Let p(x) = f(x) on the 
interval x* 2 < x < x* , where the function f(x) and the constant x* are yet to 
be determined. Using Eq. (48), we find that 

ffx 1 / 2 ) 

p(x) = 1 + ^ 1/2 ; for x* 4 < x < x* 2 , (49) 
P (x) = l + ^ + f^ for x**<x<x* 4 , (50) 

^) = 1 + i^ + i^37I + l^ for x^<x<x** , etc. (51) 

Thus we have found a family of solutions p(x) to Eq.(48), which is quite general 
since f(x) is still undetermined. 

We now show that f(x) = 1 on the interval [x* 2 ,x*). To begin, we note 
that in the stationary state, Eq. (47) implies 

In particular, for x = 1 we have 

(y + ^(l) = 2, (53) 

so that if / ~ AeP x * , then p(l) ~ 2Ae~ l3( - 1 ~ x *\ It is readily seen that p(x) ~ 
2Ae~ f3( - x ~ x *' > satisfies Eq. (52) for x* < x < 1. The same equation then leads 
to p(x) = 1 as (3 — > oo, for a;* 2 < a: < x*. We may then develop the full 
solution using Eqs. (49) to (51); normalization implies x* = 1/2. The result 
is the function p(x) plotted in Fig. 7, which is in good agreement with the 
density found via simulation of the random-neighbor model. This solution is 
discontinous at * ,..., and exhibits an integrable divergence at x — 0. 

The value of x* can be confirmed through the reasoning developed in pre- 
vious subsection for the peripheral square model, which in this case implies 
1 = 2(1 - x*) so that x* = 1/2. 



IV. NUMERICAL SIMULATION AND CRITICAL EXPONENTS 

We now compare the mean-field theory predictions with simulation results. 
We estimate the probability density p(x) on the basis of a histogram of barrier 
frequencies, dividing [0,1] into 100 subintervals. Histograms are accumulated 
after N st time steps, as required for the system (a ring of iV sites) to relax 
to the stationary state. In order to improve statistics, we average over N r 
realizations. 
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The simulation results for the original BS model are shown in Fig. 8. We 
observe qualitative agreement between nearest-neighbor (NN) and random- 
neighbor (RN) versions. (The simulation parameters are: N = 1000, N st = 
10 6 , N r = 10 3 (NN); N = 2000, N st = 10 5 , N r = 10 3 (RN).) Here and in all 
other cases, the simulation result for the random neighbor model appears to 
converge (as expected) to the mean-field prediction. In light of the discussion 
in Sec. IIIB, it is reasonable to regard the rounding of the step function as a 
finite-size effect. 

Fig. 4 presents similar results for the radical variant. We notice that the 
NN and RN versions exhibit qualitatively similar probability densities, differing 
mainly in the value of the threshold x*. (In this case we use N st = 10 6 and 
N r = 10 3 , with system sizes N = 1000 (NN) and 2000 (RN).) 

Simulation results for the centered and peripheral square versions are shown 
in Fig. 9 and 10, resp.. For the centered square, the NN and RN, probability 
densities are quite similar, differing mainly in the value of the threshold and in 
the inclination of the central portion (Fig. 9). (The simulation parameters are 
N = 2000, N st = 10 7 , N r = 500 (NN), N = 10000, N st = 10 6 , N r = 10 3 (RN).) 
On the other hand, for the peripheral square version, the nearest-neighbor and 
random- neighbor densities are somewhat different, since the latter exhibits 
various steps, while only one such step is evident in the former (see Fig. 10). 
Further study is needed to determine whether this is represents a qualitative 
difference between the two formulations, or is instead due to finite size and/or 
finite numerical resolution. (Note that for the square models we use larger 
lattices, which were necessary to observe clear singularities. The simulation 
parameters are: N = 1000, N st = 10 6 , N r = 10 3 (NN), N = 50000, N st = 10 7 , 
iV r = 20 (RN).) 

Several quantities are known to display power-law behaviour in the Bak- 
Sneppen model [1,6,18]. In particular, we studied the distribution Pj(r) that 
sucessive updated sites are separated by a distance r. In the original model, 
Pj(r) ~ r _7r , with 7r = 3.23(2) [18]. (Figures in parentheses denote uncertain- 
ties.) We performed simulations of the modified models using 10 9 time steps 
on lattices of 2000 or more sites, yielding (see Fig. 11), n = 3.22(2) for the 
original BS model; ir = 3.27(2) for the radical variant; n = 3.25(2) and 3.24(2) 
for the centered and peripheral square variants, respectively. Thus we find 
strong evidence that all of the variants introduced here belong to the same 
universality class as the original model. 

V. CONCLUSIONS 

In order to understand the implications of extremal dynamics, we propose 
several modified Bak-Sneppen models. Although different updating rules lead 
to completely different probability densities, they are always singular at one 
or more points. The step-like singularities appear in the limiting densities, 
either as (3 — > oo in the finite temperature model (on an infinite lattice), or as 
the system size N — > oo, under extremal dynamics. Thus the double limit of 
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infinite size and zero temperature is required for the BS model or it variants 
to generate a singular probability density. 

A remarkable feature common to the original model and all the variants 
considered, is that in the infinite-size limit, and under extremal dynamics, a 
certain interval D C [0, 1] has probability zero, and yet contains the extremal 
site with probability one. If we regard the density of active sites, Prob[rr e D], 
as the order parameter p, then the BS model and its variants are seen to realize 
the 'SOC limit' [4,5] p — > + . Correlations between site variables x iy Xj, ...,x n 
are zero unless one or more of these values falls in D. 

Our conclusions regarding the form of the stationary densities are based 
on mean-field analyses that are exact for the random-neighbor versions, as is 
verified numerically. We also present a pair approximation for the original 
model. An important point is that mean-field theory captures the form of the 
probability density correctly, as shown via simulations of the modified models 
on a one-dimensional lattice. The latter suggest that the critical exponents 
are independent of dynamics. 

The Bak-Sneppen model appears to be a prototype for a large universality 
class, since countless variants, beyond those presented here, are possible. We 
expect any dynamics that respects the symmetry of the original model (that is, 
spatial isotropy), and that does not introduce new conserved quantities, to have 
the same exponents as the original model. This is interesting, since the same 
critical behavior may subsist, as it were, upon stationary distributions of very 
different forms. Aside from its intrinsic interest, the question of universality is 
important for applications, since the precise form of the dynamics in a specific 
setting (e.g., evolution) is generally unknown, and probably quite different 
from that of the original model. Power laws and a singular stationary density 
only appear in the extremal dynamics limit, which may be difficult to realize 
in spatially extended natural systems. 
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FIG. 1. Original model: finite-temperature mean-field theory Eq. (6) for (3 
values as indicated ((3 inf. stands for (3 — > oo). 



1.0 



N=5 




N=20 




N=100 




N=1000 






/ 







0.0 0.2 0.4 0.6 0.8 1.0 

X 



FIG. 2. Original model on complete graph (extremal dynamics) for system sizes 
as indicated. 
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FIG. 3. Radical variant: finite-temperature mean-field theory for [3 values as 
indicated. 
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FIG. 4. Radical variant: simulation results for nearest neighbor and random 
neighbor versions compared with mean-field prediction. 
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FIG. 5. Radical variant on complete graph for system sizes as indicated. 
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FIG. 6. Centered square variant: finite-temperature mean-field theory for (3 val- 
ues as indicated. Approximation A refers to g(x) = constant, B for g{x) = ax + b, 
both at zero temperature, as explained in text. 
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FIG. 7. Peripheral square variant: finite-temperature mean-field theory for (3 
values as indicated. 
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FIG. 8. Original model: simulation results for NN and RN versions compared 
with mean-field prediction. 
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FIG. 9. Centered square variant: simulation results for NN and RN versions 
compared with mean-field prediction. 
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FIG. 10. Peripheral square variant: simulation results for NN and RN versions 
compared with mean-field prediction. 
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FIG. 11. Distribution Pj(r) of the distance r separating sucessive updated sites. 
The curves have been shifted vertically to facilitate comparison. 
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